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Abstract. We derive simple models for the dynamics of a single atom coupled to a 
cavity field mode in the absorptive bistable parameter regime by projecting the time 
evolution of the state of the system onto a suitably chosen nonlinear low-dimensional 
manifold, which is found by use of local tangent space alignment. The output field 
from the cavity is detected with a homodyne detector allowing observation of quantum 
jumps of the system between states with different average numbers of photons in the 
cavity. We find that the models, which are significantly faster to integrate numerically 
than the full stochastic master equation, largely reproduce the dynamics of the system, 
and we demonstrate that they are sufficiently accurate to facilitate feedback control of 
the state of the system based on the predictions of the models alone. 
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1. Introduction 

A two-level atom coupled to a driven and observed cavity mode exhibits a variety 
of interesting dynamics pLj, including scenarios where trajectories of the atom-cavity 
state tend to localize transiently but jump between multiple regions of phase space 
on longer timescales. Here and in related work, we loosely refer to such behaviour 
as 'bistability', and we use the term 'stable region' or 'attractor' to refer to the local 
regions of phase space, where the system spends most of its time. Bistable systems 
have potential applications as memory units and switches, and this motivates studies 
of ways to understand and control their dynamics. The phase bistable regime, where 
the system has two stable regions with different values of the phase of the cavity field, 
has been investigated in several papers [21 El HI E] , and it has been demonstrated that 
quantum jumps between the two stable regions can be observed in the photo current 
from a homodyne detector monitoring the field leaking out of the cavity [Ij. As shown in 
[6], there is also an absorptive bistable regime, for which the stable regions have different 
values of the amplitude of the cavity field mode. Quantum jump behaviour is observed 
in this case as well, but it is more complicated to obtain simple approximate descriptions 
of the dynamics due to lack of symmetry between the two stable regions [7], and we 
thus consider this regime in the following. Examples of experimental investigations of 
bistability in cavity quantum electrodynamic systems are provided in [HI El [IDl flT] . 

The time evolution of the state of a continuously monitored quantum system is 
governed by a stochastic differential equation, but it is typically a very slow process to 
integrate this equation numerically due to the large dimensionality of the Hilbert space. 
This is, in particular, a problem, if we would like to control the system dynamics through 
feedback, since, in that case, it is necessary to track the state of the system in real time. 
In many cases, it turns out that the dynamics does not explore all degrees of freedom in 
the full Hilbert space, and this opens the way to develop simple low-dimensional models, 
which can, at least approximately, predict the time evolution of the state of the system. 
One way to obtain such models is to project the system dynamics onto an affine linear 
subspace of low dimension, and this technique has turned out to be very successful in 
the case of phase bistability [5l [7] , while the results for absorptive bistability are less 
satisfactory [7]. It is, however, quite possible that improved results can be obtained by 
considering the more flexible case of projection onto a nonlinear manifold. In the present 
paper, we demonstrate that the latter approach provides simple models of absorptive 
bistable dynamics, which are sufficiently accurate to allow feedback control of the state 
of the system. 

The paper is structured as follows. In section El we introduce the system and 
integrate the stochastic master equation to provide examples of absorptive bistable 
dynamics and quantum jumps. In section El we derive reduced models of the behaviour 
of the system by first identifying a low-dimensional manifold, which captures most of 
the dynamics, and then projecting the full system dynamics onto that manifold. The 
ability of the reduced model to reproduce the results of the full model is investigated in 
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Figure 1. Two- level atom in a cavity probed with a coherent laser beam. BS is a 
beam splitter of low reflectivity, and LO is the local oscillator. 



section HI Finally, in section [5|, we demonstrate that a feedback scheme, which builds 
only on predictions of a reduced model, can be used to hold the system at one of the 
stable regions. Section [6] concludes the paper. 



2. Absorptive bistability and quantum jumps 

Our model system is a two-level atom with ground state \g) and excited state |e) coupled 
to a cavity field mode with coupling strength g as illustrated in figure [1] The cavity 
mode, which decays at a rate 2k,, is driven by a coherent laser beam, and light reflected 
from the cavity is observed with a homodyne detector. The excited state of the atom 
decays at a rate 7 by spontaneous emission, but since the emitted photons travel in 
random directions, it is so far not experimentally feasible to detect all of them with 
high efficiency. We thus assume no detection of spontaneously emitted photons in the 
following and use a density operator p to represent the state of the atom and the cavity 
field mode. 

The time evolution of p in a frame rotating with the frequency of the drive laser is 
determined by the stochastic master equation (see, for instance, [12] for a derivation) 

1 1 ^ 

dp = — — [H, p\dt + K{2apa) — a)ap — pa)a)dt + — (2o"pcr^ — aVp — pcrV)dt 

+72^ [pa^e"^ + ape-'"^ - Tr [{a^e"f' + ae"''^) p] p] dW (1) 
with Hamiltonian 

H = hAcO^a + /lAaCrV + i%o(S^cr — aa'^) + ih£{a) — a). (2) 

Here, a is the cavity field annihilation operator, a = \g){e\ is the atomic lowering 
operator, A^, is the detuning between the cavity resonance frequency and the frequency 
of the drive laser, Aa is the detuning between the atomic transition frequency and the 
frequency of the drive laser, and £ = V2kI3, where is the average number of photons 
in the probe beam arriving at the cavity input mirror per unit time and (3 is assumed 
to be real. The phase 0, which is varied experimentally by varying the relative phase of 
the probe beam and the local oscillator, determines which quadrature of the cavity field 
is detected. The x-quadrature is measured for = and the p-quadrature is measured 
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for = 7r/2. Finally, the Wiener increment dW is a Gaussian stochastic variable with 
mean and variance dt, which is related to dy, the observed homodyne photo current 
(in units of photons per time) integrated from t to t + dt divided by the square root of 
the average number of photons per unit time in the local oscillator beam, through 



In an experiment, the photo current is measured as a function of time, and we can 
eliminate dW between ([1]) and ([3]). In numerical simulations, on the other hand, we use 
a random number generator to obtain realizations of dW and integrate ([T]) directly. 

Useful insight into the dynamics predicted by ([1]) can be obtained through various 
semiclassical approximations [H [13]. We shall not pursue such models further here, 
but simply choose a set of parameters for which the dynamics of the system has been 
shown to be absorptive bistable ^j: ^c/l± = 0, Aa/7_L = 0, fi;/7_|_ = 0.1, go/'y± = \/2 
and S/'j± = 0.56, where 7_l = 7/2 is the transverse atomic decay rate. We assume 
throughout that the initial state of the system is the state with zero photons in the 
cavity and the atom in the ground state, and we use the second order derivative free 
predictor-corrector method of [H] to integrate ([1]). The expectation value of the number 
of photons in the cavity is typically below 22, and we thus truncate the basis of the 
Hilbert space of the cavity field at 59 photons, leading to a density matrix of dimension 



In figure [2] we show results for the time evolution of (a + a'^)/2 and the expectation 
value of the Pauli operators = c + a'^ and = [a^, cr]/2 for a given realization of the 
measurement noise dW. The system is seen to jump between two stable regions with 
different expectation values of the operators even though {a + ci^)/2 fluctuates more in 
the upper region than it does in the lower region. The symmetry of ([1]) dictates that 
(ay) = i{a — (T^) and — i(a — a^)/2 are both zero for the case of homodyne detection of 
the x-quadrature, while they fluctuate randomly around zero for the case of homodyne 
detection of the p-quadrature. Note that since the time average of —i{a — a^) /2 is zero for 
both of the stable regions, our ability to distinguish the regions through a measurement 
of the p-quadrature relies on the fact that the time evolution of the state is different in 
the two regions. 

3. Model reduction 

We now turn to the problem of deriving a simplified low-dimensional differential 
equation, which, ideally, contains the same dynamics as the full stochastic master 
equation ([T]). To do so, we first need to identify a suitable manifold onto which 
we can project the dynamics. Several manifold learning strategies have already been 
investigated in the literature [IHl [16] , and here we use the method of local tangent space 
alignment (LTSA) pLTj. An advantage of this method is that it optimizes the choice 
of low- dimensional space in local areas, which means that it is well suited to describe 
systems with more than one stable region. LTSA defines the manifold in terms of single 
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Figure 2. Stochastic time evolution of (a + a^)/2 (blue), (ctz) (green) and {ax) (red) 
for homodyne detection of the a;-quadrature (a) and for homodyne detection of the 
p-quadrature (b). 



points, but in order to perform the projection, we need a differentiable function, which 
relates the coordinates of all points in the low- dimensional space to the coordinates of 
the same points in the full space. This problem is solved by fitting a function to the 
points obtained from LTSA. 



3.1. Identification of the low- dimensional manifold 

In brief, the input to the LTSA algorithm is a set of vectors x*^*-*, i = 1,2, . . . , N, 
sampled with noise from an unknown d-dimensional (nonlinear) manifold embedded in 
an m-dimensional space, where m > d, and the objective is to identify the underlying d- 
dimensional manifold. For a linear manifold this is done by computing the d-dimensional 
affine subspace, which minimizes the sum of the square of the errors between the 
original vectors and the vectors projected onto the affine subspace. To tackle the more 
general case, the LTSA procedure aligns local linear structures into a global nonlinear 
manifold, where the local structures are the affine subspaces obtained by applying the 
above procedure to subsets of the A^ points. The ith subset is chosen as the k nearest 
neighbours of the ith point (including the point itself), where k is a number satisfying 
d <^ k <^ N . The output of the algorithm is the coordinates r*^*) of the points in the 
(i-dimensional space and an approximate map from the (i-dimensional space to the full 
m-dimensional space, from which it is possible to compute corrected coordinates x'-*-' of 
the points in the m-dimensional space. The map is, however, only valid in small regions 
around each r*^*\ and it is not differentiable at all points. 

In our case, we start from a set of density matrices sampled from the time evolution 
of the state of the system. Concretely, we choose the density matrix at times t = 501 7j^, 
t = 502 7j^, . . ., t = 2500 for the trajectory used to compute the results shown 
in figure [2], and we choose A; = 60 as in [7j. These density matrices are transformed 
into real column vectors x^^\ each with 14400 elements, by concatenating the real part 
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of the columns of the upper right triangular part including the diagonal followed by 
concatenation of the imaginary part of the columns of the upper right triangular part 
excluding the diagonal, i.e., 

•^n{n+l)/2 ~ ynny \^) 
m(m—l)/2+n 

Re(p«), (5) 

A''(A''+l) /2+{m— l){»Tt— 2)/2+n Mp^L), (6) 
where n = 1, 2, . . . , in (jlj) and n = 1, 2, . . . , m — 1 and m = 2, 3, . . . , in ([5]) and 
([6]). This construction method ensures that every vector x in the m-dimensional space 
corresponds to a Hermitian matrix. Furthermore, the LTSA algorithm ensures that the 
new points in the m-dimensional space are correctly normalized. There is, however, no 
guarantee that the constructed points correspond to positive semi-definite matrices, and 
here we rely on the ability of the time evolution equation to keep the state of the system 
within the physically acceptable region. Depending on the purpose of the reduced model, 
it is not necessarily optimal to minimize the projection error with respect to the dot 
product (x^'^^)'^x^^\ and one could, for instance, consider to multiply the density matrix 
elements with different weight factors [7j. We note in particular that {x^'^'^)'^x^^^ is equal 
to Tr(p(*)p(-^^) if a factor of \/2 is included on the right hand side of ([5]) and ([6]), but 
we have avoided to do so in the following, because we obtain better results without the 
factor -\/2. 

Having obtained a set of points t^^^ in the low- dimensional space and the 
corresponding coordinates x'-*-' in the full space, we next construct a map from the 
low-dimensional space to the full space via fitting. We need to compute one fit for each 
of the m = 14400 coordinates in the full space, and to make this procedure practical, 
we use the same fitting model for all the coordinates and choose this model to be linear 
in the fitting parameters, i.e., we assume a map of form x = cf{T), where x is a vector 
in the m-dimensional space, c is an m x r matrix of fitting parameters, and / is an r x 1 
vector, whose elements fj are arbitrary functions of the coordinates r = (ri, T2, . . . , r^)^ 
in the (i-dimensional space. To ensure that x is correctly normalized for all r, we 
choose /i(t) = 1 and minimize — cf (t^^^))^ (x^^^ ~c/(r'^*^)) under the constraint 

v'^c = (1, 0, . . . , 0), where f is an m x 1 vector, whose ith entry is one ii i = n{n + l)/2 
for some n G {1, 2, . . . , A^} and zero otherwise (i.e., v'^x = Tr(p), where p is the density 
matrix corresponding to the vector x). The result is 

c = c+-^v[{l,0,...,0)-v^c], (7) 

where c = [{z'^z)~'^z'^y]'^ , with Zij = /j(T^*)) and yij = x^j\ is the standard linear least 
squares result without constraints. 
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3.2. Projection of the dynamics onto the identified manifold 
The result of the last subsection is a relation of form 

Pi.'^) = ^Cjfj{Ti,T2,...,Td), (8) 

i 

where Cj is the matrix obtained by applying the inverse of (HUH]) to the jth column of c. 
To project the stochastic master equation onto the manifold defined by ([8]), we follow 
the derivation in |7j . We would like to interpret dp as a vector, but dp only transforms as 
a vector if (diy)^ = 0, i.e., if we use Stratonovich calculus, and we thus rewrite ([1]) into 
Stratonovich form dp = A[p]dt + i?[p] odW^. Starting from a point p(r) on the manifold, 
we then project dp(r) onto the tangent space of the manifold at that point, which is 
spanned by the d vectors dp{T)/dTi, using the dot product {pa, Pb) = T^t^{paPb), i-e., 

+ EE(«-%TJflW^)]^}^odH^ (9) 

i j ^ J J 'I- 

where g = gir) is the metric tensor with elements 

dp dp\ _^^dfp dfg 



Combining this relation with 

dpW = E^°d-.. (11) 

i 

we obtain an expression for the time evolution of 

B[p{r)]^\ o dW. (13) 



+ 5^(r%Tr| 



3 

To simplify the notation, we define vectors vi and f 2 with elements 

{vi)j = v^Tr[(ate*'^ + de-''^)cj], (14) 

{V2)j = 2KRe{Tr[(a^e^'^ + de-''^)de-'^cj]}, (15) 

and matrices Mg, Mh and Mi with elements 

{Mg),, = Tr(Q9), (16) 

{MH)ij = 2Re[Tr(Mc,Q)] + 27Tr(a_c,a+Ci), (17) 

{Mi),j = V2^Tr[(cja^e^'^ - ae-^'^cj)ci], (18) 

where M = —iH/h — K.a)d — Kci^e"^*'^ — 70"+ (T_. Finally, we insert ^[p(t)] and i?[p(r)] 
into (JT2I) and convert back to ito form to obtain 

dn = ai [t] dt + hi [t] dW, (19) 
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where 



and 



\dTkdTj dTj dTk dTk dTj orkdrj ^ / 

4 E + f M,^) ..MM^] (20) 



h[r] = Y.(g-% {^-£-M,f - vlf^MJ^ . 



(21) 



Integrating the low-dimensional equation ([T^ . we can now approximately predict the 
time evolution of p through ([8]). 



4. Performance of the reduced models 



Since we would like to use the reduced models to predict the state of the system in a 
feedback scheme, we should check the performance of the reduced models by generating a 
realistic photo current using the full stochastic master equation and then use that photo 
current to integrate the reduced models. Examples of this procedure, using polynomials 
as fitting models, are provided in figureO We have plotted (a + a)) /2, because this is the 
quantity we need to predict in the feedback scheme proposed in the next section. For 
the case of homodyne detection of the p-quadrature we have used a weighted linear least 
squares method to compute c in ([7]) to reduce the effect of outliers. The precise initial 
state of the reduced models is not important, because the observed value of the photo 
current quickly drags the reduced model to the correct state, and we have thus chosen 
r = (0, 0, ... , 0)"^ (the average of r^*^) for simplicity. In case of instability, we have reset 
the reduced model to r = (0, 0, ... , 0)^ whenever the program returns a non-determined 
value of r. 

The figure shows that simple low-dimensional models are able to provide accurate 
predictions for the value of (a + a)) /2 for the case of homodyne detection of the x- 
quadrature. For homodyne detection of the p-quadrature, the two-dimensional model 
obtained by using a second order polynomial as fitting model is observed to largely 
reproduce the time evolution of (a -|- a^)/2, but the details differ. In particular, the 
reduced model does not reproduce the highest values of (a -|- a^)/2. To improve the 
agreement between the full and the reduced model, one could try to increase the order of 
the fitting polynomial or use a higher- dimensional model. For a fourth order polynomial 
in two dimensions, the reduced model is able to reach the highest values of (a + a^)/2 
and the value of (a + a^)/2 in the lower state is also more accurate, but the model 
predicts false jumps to the lower state, which is undesirable in a feedback scheme. False 
jumps are not observed in the figure for the four- dimensional model with a second order 
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Figure 3. Time evolution of (a + at)/2 obtained from the full stochastic master 
equation (blue) and from various reduced models (green) assuming the same photo 
current. For (a) and (d) the fitting model is a second order polynomial in two 
dimensions, i.e., / — (1, ri, T2, t^, Tir2, t|)"^, for (b) and (e) it is a fourth order 
polynomial in two dimensions, and for (c) and (f) it is a second order polynomial 
in four dimensions. The left figures are for homodyne detection of the x-quadrature, 
and the right figures are for homodyne detection of the p-quadrature. 

polynomial as fitting model, and we thus use this latter model to predict the state of 
the system in the next section. 

The reason why we obtain very accurate results for homodyne detection of the x- 
quadrature is that we use the actual photo current dy obtained from the full stochastic 
master equation to drive the reduced model. This means that the back action of the 
measurement on the system is large whenever the predicted value of the measured 
quadrature differs significantly from the value obtained from the full stochastic master 
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equation. It is, in fact, a much harder test of the performance of the reduced models 
to check whether they are able to predict the value of quantities that are not related in 
a simple way to the observed quadrature such as the value of (d + d^)/2 for homodyne 
detection of the p-quadrature. The above results thus indicate that the low-dimensional 
models actually capture most of the full dynamics of the system. As expected, we also 
find that the reduced models for homodyne detection of the p-quadrature provide more 
accurate results for —i{a — a^) /2 than for {a + a^) /2. 

5. Stabilization of one attractor through feedback control 

A natural feedback scheme to hold the system within one of the stable regions is to 
increase or decrease the intensity of the drive laser depending on whether the value of 
(d + d^)/2 is below or above some suitably chosen target value xq. This is achieved by 
adding a proportional feedback term 

dpf|3 = Spe(t)[d - d^p]7_Ldt, e(t) = (d + d"'')/2 - a;o, (22) 

to the stochastic master equation ([T]), where Sp is a parameter determining the strength 
of the feedback. Integration of the resulting equation confirms that the feedback term 
has the desired effect, and it is possible to decrease the standard deviation of e{t) to 
a value, which is very small compared to typical fluctuations of (d + d^)/2 without 
feedback. The question is now whether this is still the case if we use a reduced model 
(with the feedback term included) to predict the value of e{t). 

Example trajectories are shown in figure |H Comparing these trajectories to those 
in figure [HI it is apparent that the feedback term affects the time evolution, and that 
the system stays close to the upper or the lower stable region for the considered values 
of Xq. For the upper stable region, (d + d^')/2 is seen to oscillate with a relatively large 
amplitude, which reflects the fact that the upper stable region is relatively broad as 
observed in figure [2l For the lower stable region, (d + d^")/2 is roughly constant except 
for sudden spikes. By choosing a value of xq, which is slightly below the average value 
of (d + d"'")/2 predicted by the reduced model, we can ensure that the feedback term 
almost always acts to decrease (d + d)) /2. This keeps the full model away from the 
transition region, and as seen in the figure the spikes tend to point towards negative 
values of (d + d^)/2 rather than towards larger positive values of (d + d''')/2 as observed 
for higher values of xq. 

For the case of homodyne detection of the p-quadrature, we note that the 
predictions of the reduced model for xq = 3.5 fluctuate less than the results obtained 
from the full stochastic master equation, while the fluctuations of the full and the 
reduced model are approximately the same for the case of homodyne detection of the 
x-quadrature. This is because we use the reduced model to evaluate the error e(t), 
which means that the feedback term always acts to reduce e{t) for the low-dimensional 
model. For the full model, on the other hand, e{t) may have the opposite sign, in which 
case the full model is pushed away from xq. The resulting change in the photo current 
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Figure 4. Time evolution of (a + dt)/2 obtained from the full stochastic master 
equation (blue) for homodyne detection of the a:-quadrature (a) and homodyne 
detection of the p-quadrature (b) for the same noise realization as in figure [31 but 
with the feedback term in (|22p included. The error e{t) is computed from the value 
of (a + a^)/2 predicted by the four-dimensional reduced model with a second order 
polynomial as fitting model (green), and as in the last section we have used the 
photo current obtained from the full stochastic master equation to integrate the low- 
dimensional model. In (a), Xq = 3.8 and Sp = 0.75 for the upper curves and xq = 0.4 
and Sp — 1 for the lower curves. In (b), xq = 3.5 and Sp = 0.75 for the upper curves 
and xq = 0.5 and Sp — I for the lower curves. 



drives the reduced model in the same direction as the full model, but this mechanism is 
more efficient in the case of homodyne detection of the x-quadrature than for homodyne 
detection of the p-quadrature. 

Discrepancies between the full model and the reduced model lead to a feedback of 
noise into the system, and even though a large value of Sp may reduce the variance of e{t) 
for the low- dimensional model, we observe that the predictions for (a + a)) /2 obtained 
from the full model fluctuate over a range that is broader than the distance between 
the upper and the lower stable region if Sp is chosen too large. The power spectrum 
of the time evolution of (a -|- d))/2 for a trajectory computed from the full stochastic 
master equation with homodyne detection of the p-quadrature and the power spectrum 
of the difference between the predictions of the reduced model and the full model in 
figure [5]^a) show that the relative error of the reduced model in predicting (a -|- d)) /2 is 
smaller at low frequencies. This appears because the reduced model is able to predict 
quantum jumps of the system but does not capture the details of the dynamics in the 
stable regions, and it suggests that it might be an advantage to mainly feed back the low 
frequency behaviour, which can be achieved by adding an integral term to the controller 



dPfb 



Spe{t) + Si I exp {—({t — t')) e{t')dt' 



a^p]7xdt, (23) 



where Si and ( are constants. 

Rough estimates of reasonable choices of Sp, Si and ( can be obtained as follows. 
If the measurement noise is turned off by setting dW = in the full stochastic master 
equation, the state of the system decays to steady state, which is an incoherent mixture 
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Figure 5. (a) Power spectrum of a trajectory corresponding to the one in figure [Sj^f), 
but integrated to tj± — 4000. The blue curve is the power spectrum of (a + a^)/2 
obtained from the full stochastic master equation, the green curve is the power 
spectrum of the difference between the reduced and the full model, and uj is the angular 
frequency. In both cases we have subtracted the mean value before computing the 
power spectrum, (b) Plots of |wJ(C — iu!)Q\ (dashed) and of the norm (solid) and the 
phase (dotted) of the loop transfer function K{iuj). 



of the states corresponding to the upper and the lower stable regions. (Note that dW 
in equation (fT9l) may be nonzero, since we still require that the the photo current is the 
same for the reduced and the full model.) The behaviour of the system when a small 
input drive field term dp„ = u(t)[d — d'^ , p]'y±dt is added can then be investigated by 
linearizing ([T]) and (fT9|) around the steady state point, which leads to an equation of 
form 



dt 



5r 
6p 



C 



5r 
6p 



+ Qu{t), 



(24) 



where 6t is the deviation of r from the steady state value, 6p is a vector of the deviations 
of the density matrix elements from their steady state values with one diagonal element 
omitted, C is a 14403 x 14403 matrix and Q is a 14403 x 1 vector. Choosing xq to be 
the steady state value of (a + a^)/2 obtained from the reduced model, we can write the 
error as 

^ 6t 



e{t) 



6p 



(25) 



where is a 14403 x 1 vector for which all but the first four elements are zero. Integrating 
the linear equation (^^, we have an expression for the time evolution of e(t), which 
we insert into (1231) to obtain the feedback term for a given input u{t). In a closed 
loop setting, the feedback term is used as input, and the system behaviour is thus 
characterized by the loop transfer function 



K{s) = - 



Sp + 



C + s 



(26) 



which is the Laplace transform of the coefficient in square brackets in (1231) divided 
by the Laplace transform of u{t). For s = iu, the norm of K is the loop gain at 
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Figure 6. Same as in figure [H but for the feedback term in ([^5]) . The parameters are 
Sp = 0, Si — 0.51 and ( = 2tt x 0.1 j± for the case of homodyne detection of the 
x-quadrature and Sp — 0, Si — 0.15 7^ and C, = 2tt x 0.02 7^ for tlic case of homodyne 
detection of the p-quadrature. 

angular frequency u, and according to the power spectra in figure [5]^a) this should be 
close to zero for angular frequencies above approximately 2tt x 0.02 7^ = 0.126 7^. 
The norm of vJ(C — iu!)~^Q is determined completely by the system and is plotted in 
figure [5](b) for the case of homodyne detection of the p-quadrature. Since this factor is 
substantially different from zero for a range of angular frequencies above 2tt x 0.02 7^, 
we choose Sp = 0. To ensure that \si{( + iuj)~^\ has a large negative derivative for 
u; ~ 27r X 0.02 7_l, we set the angular cross-over frequency to = 27r x 0.02 7_l. Finally, 
we choose Si = V2C{\v'^{C - iCy^Ql)-^ = 0.26 7^ such that \K{iC) \ = 1. The resulting 
norm and phase of the loop transfer function are also plotted in figure [5]^b). The phase 
is seen to be well above — tt at the cross-over frequency, and we thus expect the feedback 
to be stable. One should, however, not read too much into (126|) as the above derivation 
is very crude. A similar analysis for the case of homodyne detection of the x-quadrature 
leads to the parameters Sp = 0, Sj = 0.51 7j_ and C = 27r x 0.1 7j_ and suggests that the 
feedback may be unstable for Sj > 0.7 7_l. 

Keeping Sp and ( fixed and searching in the neighbourhood of the above values of 
Si, we find that Si ~ 0.51 7^ is close to optimal for the case of homodyne detection of the 
x-quadrature. A too small value of Si (for instance, below 0.25 7^) leads to increased 
fluctuations in the predictions of the full model, because the feedback is insufficient 
to keep the system within the stable region, and a too large value of Sj (for instance, 
above 1 7^) tends to course instability. For the case of homodyne detection of the p- 
quadrature, we obtain improved results by decreasing Si to 0.15 7^. Trajectories for these 
parameters are shown in figure [61 For the upper stable region the standard deviation of 
the difference between the reduced and the full model relative to the standard deviation 
of the trajectory obtained from the full model for the time interval from t = 100 7j^ 
to t = 750 7j^ is reduced relative to the value obtained for the trajectories in figure HI 
For homodyne detection of the x-quadrature it decreases from 0.71 to 0.46, and for 
homodyne detection of the p-quadrature it decreases from 1.01 to 0.82. This confirms 
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that the integral term feeds back less noise. On the other hand, the integral term is 
less efficient in keeping {a + a))/2 close to the desired value, and the standard deviation 
of (a + a^)/2 obtained from the full model is observed to increase. The conclusion is 
the same for the lower stable region for homodyne detection of the x-quadrature, but 
for homodyne detection of the p-quadrature we observe a decrease in both the standard 
deviation of the difference between the reduced and the full model and in the standard 
deviation of {a + a)) /2 obtained from the full model. 

6. Conclusion 

In conclusion, we have shown that it is possible to derive simple low-dimensional models 
for the dynamics of a single atom interacting with a cavity field mode in the absorptive 
bistable regime, and we have demonstrated that the models can be used to construct a 
feedback scheme, which is able to hold the system within one of the stable regions. This 
is an important result, because it is unrealistic to integrate the full stochastic master 
equation in real time. 

The suggested feedback scheme relies on predictions of the expectation value of 
the ,x-quadrature of the cavity field, and wc have considered both the case of homodyne 
detection of the x-quadrature of the output field from the cavity and homodyne detection 
of the p-quadrature. The former case is relatively easy to handle, because the estimated 
quantity is directly related to the observed quantity. For homodyne detection of the 
]9-quadrature, on the other hand, the expectation value of the x-quadrature has to be 
inferred from the precise time evolution of the p-quadrature, and the reduced model has 
to do significantly more work. It is thus promising for the method that we also obtain 
reasonable results in this case. 

There are many degrees of freedom in the modelling procedure and, in general, it 
may require some trial and error to find reduced models that are able to reproduce the 
system dynamics with sufficient accuracy. In a feedback scheme, one should concentrate 
on optimizing the ability of the reduced model to predict the quantity that determines 
the feedback, since errors in this quantity lead to a feedback of noise into the system, 
which limits the performance of the feedback. 

A further line of research could be to find systematic methods to optimize the 
models. One could, for instance, consider different ways to construct the vectors used 
as input to the local tangent space alignment procedure, which corresponds to different 
criteria for the optimal choice of low-dimensional manifold. It is also possible that 
improved models could be obtained by choosing other kinds of fitting models than 
polynomials of low order. The selection of the density operators used to compute the 
low-dimensional manifold could be adjusted according to the purpose of the model. In 
the case of application of feedback, one could, for instance, include more points at one 
attractor than the other in order to obtain a better description of the dynamics in the 
neighbourhood of the attractor we intend to stabilize. 
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